fig = plt.figure()
plt.plot(np.log(np.abs(z_left)),np.log(rho_left), "r-", label = "$z \\in [-L/2,0]$")
plt.plot(np.log(np.abs(z_right)),np.log(rho_right), 'b-', label = "$z \\in [0,L/2]$")
plt.legend()
plt.xlabel("$log|z|$")
plt.ylabel("$log|\\rho|$")
plt.show()
plt.plot(np.log(-z_left[::-1]),np.log(rho_avgd))
plt.legend()
plt.xlabel("$log|z|$")
plt.ylabel("$log|\\rho|$")
plt.show()
#ADDITIONAL:
#Calculate potential
phi = GF.fourier_potentialV2(rho,L)
#Calculate Acceleration Field on Mesh:
a_grid = NB.acceleration(phi,L)
--------------------------------------------------------------------------- NameError Traceback (most recent call last) /home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb Cell 7' in <cell line: 2>() <a href='vscode-notebook-cell:/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb#ch0000006?line=0'>1</a> fig = plt.figure() ----> <a href='vscode-notebook-cell:/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb#ch0000006?line=1'>2</a> plt.plot(np.log(np.abs(z_left)),np.log(rho_left), "r-", label = "$z \\in [-L/2,0]$") <a href='vscode-notebook-cell:/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb#ch0000006?line=2'>3</a> plt.plot(np.log(np.abs(z_right)),np.log(rho_right), 'b-', label = "$z \\in [0,L/2]$") <a href='vscode-notebook-cell:/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb#ch0000006?line=3'>4</a> plt.legend() NameError: name 'z_left' is not defined
<Figure size 432x288 with 0 Axes>
#hist, xedges, yedges = np.histogram2d(stars_x,stars_v,bins = [50,50])
#plt.imshow(hist,extent = (x_min,x_max,v_min,v_max), cmap = cm.hot, aspect = (x_max-x_min)/(v_max-v_min))
#plt.show()
#Re-center the system
z = np.linspace(-L/2,L/2,N)
max_index = 620
print(z[max_index])
stars_x_new = stars_x - z[max_index] #centroid_z
nbins = 100
hist, xedges, yedges, image = plt.hist2d(stars_x_new, stars_v,
bins = [nbins,nbins],
range = [[-L/2, L/2],[np.min(stars_v),np.max(stars_v)]],
#cmax = 250,
cmap = cm.hot)
plt.colorbar()
plt.show()
# print(hist)
# for i in range(len(hist)):
# for j in range(len(hist)):
# hist[i,j] = int(hist[i,j])
# print(hist)
def Abel(x_range: tuple, v_range: tuple,f):
# f is a 2d-array. (A 2d histogram)
x_min,x_max = x_range
#v_min,v_max = v_range
n_rows,n_cols = np.shape(f)
#v_s = np.linspace(v_min,v_max,n_rows)
x_s = np.linspace(x_min,x_max,n_cols)
#dv = v_s[1]-v_s[0]
dx = x_s[1]-x_s[0]
holder = np.array([])
for i in range(n_rows):
#v = v_s[i]
sum = 0
for j in range(n_cols):
#x = dx*j #x_s[j]. Note we only add up from x=0 to x->inf
sum+= f[i][j]*dx
# r = np.sqrt(v**2+x**2)
# dr = np.sqrt(dv**2+dx**2)
#
# term1 = f[i][j]*r
# term2 = np.sqrt(r**2-v**2)
# sum+=dr*term1/term2
holder = np.append(holder,sum)#2*sum)
return holder
np.savetxt("hist.csv",hist,fmt ='%i',delimiter = ",")
hist = np.loadtxt("hist.csv",dtype = int,delimiter = ',')
v_min,v_max = np.min(stars_v),np.max(stars_v)
abel_transform = Abel((-L/2,L/2),(v_min,v_max),hist)
#print(abel_transform)
v_array = np.linspace(v_min,v_max,len(abel_transform))
plt.plot(v_array,abel_transform,'b--',marker = '.')
plt.show()
0.24124124124124124
import Analysis
#[ ..., [r,m,Num_bosons,sigma,Num_stars],...]
Args = [
[0.5,1.0,0,1,10000],
[0.5,1.0,10000,1,10000],
[1,0.5,20000,1,10000],
[5,0.1,100000,1,10000],
[10,0.05,200000,1,10000],
[50,0.01,1000000,1,10000],
[0.5,1.0,10000,1,0]
]
for args in Args:
print("----------New Analysis--------")
print(
f"r = {args[0]}",
f"mu = {args[1]}",
f"Num_bosons = {args[2]}",
f"sigma = {args[3]}",
f"Num_stars = {args[4]}"
)
Analysis.analysis(*args)
/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis ----------New Analysis-------- r = 0.5 mu = 1.0 Num_bosons = 0 sigma = 1 Num_stars = 10000
/home/boris/Documents/Research/Coding/OneD/WaveNonDim.py:128: RuntimeWarning: invalid value encountered in true_divide F_s = F_s/Norm_const
0.002002002002002068
chi^2 = 1069130.1514376372
chi^2 = 331220.8341718659
----------New Analysis-------- r = 0.5 mu = 1.0 Num_bosons = 10000 sigma = 1 Num_stars = 10000
0.002002002002002068
chi^2 = 1923005.6906601514
chi^2 = 1931769.6458091086
chi^2 = 500753.0354785936
----------New Analysis-------- r = 1 mu = 0.5 Num_bosons = 20000 sigma = 1 Num_stars = 10000
0.002002002002002068
chi^2 = 2025205.074171987
chi^2 = 1022913.9619586606
chi^2 = 617435.4881954561
----------New Analysis-------- r = 5 mu = 0.1 Num_bosons = 100000 sigma = 1 Num_stars = 10000
0.002002002002002068
chi^2 = 916046.8264148047
chi^2 = 851932.3919991729
chi^2 = 328193.44588331267
chi^2 = 330860.4024658548
----------New Analysis-------- r = 10 mu = 0.05 Num_bosons = 200000 sigma = 1 Num_stars = 10000
0.0
chi^2 = 1738057.7770754427
chi^2 = 1716889.4310029696
chi^2 = 284629.91746737063
chi^2 = 291832.4249495785
----------New Analysis-------- r = 50 mu = 0.01 Num_bosons = 1000000 sigma = 1 Num_stars = 10000
0.002002002002002068
chi^2 = 309143.8474700038
chi^2 = 321876.3329614536
chi^2 = 494952.2346602648
/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/Analysis.py:437: RuntimeWarning: invalid value encountered in power og = a0/((z**a2) * (z+a1)**a3)
chi^2 = 367295.279247727
----------New Analysis-------- r = 0.5 mu = 1.0 Num_bosons = 10000 sigma = 1 Num_stars = 0
import numpy as np
import matplotlib.pyplot as plt
My_Package_PATH = "/home/boris/Documents/Research/Coding"
import sys
sys.path.insert(1, My_Package_PATH)
import OneD.NBody as NB
z = np.linspace(-1,1)
x = 0
v = 1
star = NB.star(0,1,x,v)
dt = 0.1
t = 0
i = 0
while t < 2:
plt.plot(star.x,0,'ro')
plt.xlim(-1,1)
plt.show()
star.x -= v*dt
star.reposition(2)
t += dt
div = -1.0, remainder = 0.9 new z = -0.09999999999999998
div = -1.0, remainder = 0.8 new z = -0.19999999999999996
div = -1.0, remainder = 0.7000000000000001 new z = -0.29999999999999993
div = -1.0, remainder = 0.6000000000000001 new z = -0.3999999999999999
div = -1.0, remainder = 0.5000000000000001 new z = -0.4999999999999999
div = -1.0, remainder = 0.40000000000000013 new z = -0.5999999999999999
div = -1.0, remainder = 0.30000000000000016 new z = -0.6999999999999998
div = -1.0, remainder = 0.20000000000000018 new z = -0.7999999999999998
div = -1.0, remainder = 0.1000000000000002 new z = -0.8999999999999998
div = -1.0, remainder = 2.220446049250313e-16 new z = -0.9999999999999998
div = -2.0, remainder = 0.9000000000000001 new z = 0.9000000000000001
div = 0.0, remainder = 0.8000000000000002 new z = 0.8000000000000002
div = 0.0, remainder = 0.7000000000000002 new z = 0.7000000000000002
div = 0.0, remainder = 0.6000000000000002 new z = 0.6000000000000002
div = 0.0, remainder = 0.5000000000000002 new z = 0.5000000000000002
div = 0.0, remainder = 0.40000000000000024 new z = 0.40000000000000024
div = 0.0, remainder = 0.30000000000000027 new z = 0.30000000000000027
div = 0.0, remainder = 0.20000000000000026 new z = 0.20000000000000026
div = 0.0, remainder = 0.10000000000000026 new z = 0.10000000000000026
div = 0.0, remainder = 2.498001805406602e-16 new z = 2.498001805406602e-16